function [x_nodes] = ex33_mesh_generation2(interval_number,r_int,L)
%function [x_nodes] = ex33_mesh_generation2(interval_number,r_int,L)
% x_nodes:         nodal coordinates 
% interval_number: number of intervals (it must be a multiple of 3)
% r_int:           internal radius
% L:               dimension of the square

if mod(interval_number,3) ~= 0
    disp('Error: the number of intervals must be a multiple of 3');
end


i_n     = interval_number /3;
r_lim = (L - r_int)/3 + r_int;


% mesh1
k       = 1;
dtheta  =  pi/2/(i_n*2);
h       = (L - r_int)/3/i_n;
x_mesh1 = [];
for i=1:(i_n*2+1)
   theta = dtheta * (i-1);
   j=i_n+1;
   r = r_int +(j-1)*h;
   x = r * cos(theta);
   y = r * sin(theta);
   x_mesh1(k,:) = [x y];
   k = k+1;
end


N = i_n+1;
l = (L - r_int)/3;
a = (l-3/4*h-h/4*N)/(N*N-1);
b = h/4-a;

aux = 1;

for j=0:N
    r = r_int +a*(j)*(j)+b*(j);
    dtheta = pi/2/(aux*(N+1-j)+i_n*2);
    for i=1:(aux*(N+1-j)+i_n*2+1)
        theta = dtheta * (i-1);
        x = r * cos(theta);
        y = r * sin(theta);
        x_mesh1(k,:) = [x y];
        k = k+1;
    end
end





% mesh2
k       = 1;
dtheta  =  pi/2/(i_n*2);
x_mesh2 = [];
for j=1:(i_n+1)
    theta = dtheta * (j-1);
    x1    = r_lim * cos(theta);
    y     = r_lim * sin(theta);
    hx    = (L-x1) / (i_n*2);
    for i=1:(i_n*2)
        x = x1 + hx*(i);
        x_mesh2(k,:) = [x y];
        k = k+1;
    end    
end

% mesh3
k       = 1;
dtheta  =  pi/2/(i_n*2);
x_mesh3 = [];
for j=1:(i_n+1)
    theta = pi/4 + dtheta * (j-1);
    y1    = r_lim * sin(theta);
    x     = r_lim * cos(theta);
    hy    = (L-y1) / (i_n*2);
    for i=1:(i_n*2)
        y = y1 + hy*(i);
        x_mesh3(k,:) = [x y];
        k = k+1;
    end    
end


% mesh4
x1      = r_lim * cos(pi/4);
y1      = r_lim * sin(pi/4);
hx      = (L - x1) / (i_n*2);
hy      = (L - y1) / (i_n*2);
k       = 1;
x_mesh4 = [];
for i=1:(i_n*2)
    x = x1 + hx *(i);
    for j=1:(i_n*2)
        y = y1 + hy *(j);
        x_mesh4(k,:) = [x y];
        k = k+1;
    end
end


x_nodes = [ x_mesh1
            x_mesh2;
            x_mesh3;
            x_mesh4];
            
x_nodes = unique(x_nodes,'rows');


